tic
clear;
load 'C:\Users\scott\Dropbox\Spatial Probit PSRM\PSRM_R&R\App_Code\Results\hb_ksg_RR_rho_phi.mat'

%Immediate Effects 

unit_1 = 1021;%Senegal in 1993
unit_2 = 1062;%1992
unit_3 = 1103;
unit_4 = 1145;
unit_5 = 1187;
unit_6 = 1229;
unit_7 = 1271; 
    
X = XNT;
        

mult_fx = inv(i_nt - pt(10,1)*TL - pt(11,1)*WNT);
SIGMA = mult_fx*mult_fx';
rf_cut_full = mult_fx*(pt(1,1) + pt(2,1)*X(:,1)+ pt(3,1)*X(:,2)+ pt(4,1)*X(:,3)+ pt(5,1)*X(:,4)+ pt(6,1)*X(:,5)+ pt(7,1)*X(:,6)+ pt(8,1)*X(:,7)+ pt(9,1)*X(:,8)); 


%1st period 

p_1 = normcdf(rf_cut_full(unit_1),0,sqrt(SIGMA(unit_1,unit_1)));

p1_1 = mvncdf([-Inf,-Inf], [rf_cut_full([unit_1]), rf_cut_full(unit_2)],0,SIGMA([unit_1, unit_2],[unit_1,unit_2])); 

    cp_1 = p1_1/p_1;

p_0 = normcdf(-rf_cut_full(unit_1),0,sqrt(SIGMA(unit_1,unit_1)));

p1_0 = mvncdf([rf_cut_full([unit_1]),-Inf], [Inf, rf_cut_full(unit_2)],0,SIGMA([unit_1, unit_2],[unit_1, unit_2])); 

    cp_0 = p1_0/p_0;

cp_diff(1,1) = cp_1-cp_0;

%2nd period 

p_11 = p1_1; 
p1_11 = mvncdf([-Inf,-Inf,-Inf], [rf_cut_full([unit_1]), rf_cut_full(unit_2), rf_cut_full(unit_3)],0,SIGMA([unit_1, unit_2, unit_3],[unit_1,unit_2,unit_3])); 

    cp_11 = p1_11/p_11;

p_10 = mvncdf([-Inf, rf_cut_full(unit_2)], [rf_cut_full(unit_1), Inf],0,SIGMA([unit_1, unit_2],[unit_1,unit_2])); 
p1_10 = mvncdf([-Inf,rf_cut_full(unit_2),-Inf], [rf_cut_full([unit_1]), Inf, rf_cut_full(unit_3)],0,SIGMA([unit_1, unit_2, unit_3],[unit_1,unit_2,unit_3]));

    cp_10 = p1_10/p_10;
   
p_00 = mvncdf([rf_cut_full([unit_1]),rf_cut_full(unit_2)], [Inf, Inf],0,SIGMA([unit_1, unit_2],[unit_1,unit_2])); 
p1_00 = mvncdf([rf_cut_full(unit_1),rf_cut_full(unit_2),-Inf], [Inf, Inf, rf_cut_full(unit_3)],0,SIGMA([unit_1, unit_2, unit_3],[unit_1,unit_2,unit_3]));

    cp_00 = p1_00/p_00;
    
cp_diff(1,2) = cp_11-cp_00;
cp_diff(2,2) = cp_10-cp_00;



%3rd period 

p_111 = p1_11; 
p1_111 = mvncdf([-Inf,-Inf,-Inf,-Inf], [rf_cut_full([unit_1]), rf_cut_full(unit_2), rf_cut_full(unit_3),rf_cut_full(unit_4)],0,SIGMA([unit_1, unit_2, unit_3, unit_4],[unit_1,unit_2,unit_3, unit_4])); 

    cp_111 = p1_111/p_111;
    
p_110 = mvncdf([-Inf, -Inf, rf_cut_full(unit_3)], [rf_cut_full(unit_1),rf_cut_full(unit_2), Inf],0,SIGMA([unit_1, unit_2, unit_3],[unit_1,unit_2, unit_3])); 
p1_110 =  mvncdf([-Inf,-Inf,rf_cut_full(unit_3),-Inf], [rf_cut_full([unit_1]), rf_cut_full(unit_2), Inf ,rf_cut_full(unit_4)],0,SIGMA([unit_1, unit_2, unit_3, unit_4],[unit_1,unit_2,unit_3, unit_4])); 

    cp_110 = p1_110/p_110;
 
p_100 = mvncdf([-Inf, rf_cut_full(unit_2), rf_cut_full(unit_3)], [rf_cut_full(unit_1), Inf, Inf],0,SIGMA([unit_1, unit_2, unit_3],[unit_1,unit_2, unit_3])); 
p1_100 = mvncdf([-Inf,rf_cut_full(unit_2),rf_cut_full(unit_3),-Inf], [rf_cut_full([unit_1]), Inf, Inf ,rf_cut_full(unit_4)],0,SIGMA([unit_1, unit_2, unit_3, unit_4],[unit_1,unit_2,unit_3, unit_4])); 

    cp_100 = p1_100/p_100;

p_000 = mvncdf([rf_cut_full(unit_1), rf_cut_full(unit_2), rf_cut_full(unit_3)], [Inf, Inf, Inf],0,SIGMA([unit_1, unit_2, unit_3],[unit_1,unit_2, unit_3]));
p1_000 = mvncdf([-Inf,rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4)], [rf_cut_full([unit_1]), Inf, Inf ,Inf],0,SIGMA([unit_1, unit_2, unit_3, unit_4],[unit_1,unit_2,unit_3, unit_4])); 

    cp_000 = p1_000/p_000;

cp_diff(1,3) = cp_111-cp_000;
cp_diff(2,3) = cp_110-cp_000;
cp_diff(3,3) = cp_100-cp_000;

%4th period 

p_1111 = p1_111; 
p1_1111 = mvncdf([-Inf,-Inf,-Inf,-Inf,-Inf], [rf_cut_full([unit_1]), rf_cut_full(unit_2), rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5) ],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5],[unit_1,unit_2,unit_3, unit_4, unit_5])); 

    cp_1111 = p1_1111/p_1111;
    
p_1110 = mvncdf([-Inf,-Inf,-Inf,rf_cut_full(unit_4)], [rf_cut_full([unit_1]), rf_cut_full(unit_2), rf_cut_full(unit_3), Inf ],0,SIGMA([unit_1, unit_2, unit_3, unit_4],[unit_1,unit_2,unit_3, unit_4])); 
p1_1110 = mvncdf([-Inf,-Inf,-Inf,rf_cut_full(unit_4),-Inf], [rf_cut_full([unit_1]), rf_cut_full(unit_2), rf_cut_full(unit_3), Inf, rf_cut_full(unit_5)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5],[unit_1,unit_2,unit_3, unit_4, unit_5])); 

    cp_1110 = p1_1110/p_1110;
    
p_1100 = mvncdf([-Inf,-Inf, rf_cut_full(unit_3),rf_cut_full(unit_4)], [rf_cut_full([unit_1]), rf_cut_full(unit_2), Inf, Inf ],0,SIGMA([unit_1, unit_2, unit_3, unit_4],[unit_1,unit_2,unit_3, unit_4])); 
p1_1100 = mvncdf([-Inf,-Inf,rf_cut_full(unit_3),rf_cut_full(unit_4),-Inf], [rf_cut_full([unit_1]), rf_cut_full(unit_2), Inf, Inf, rf_cut_full(unit_5)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5],[unit_1,unit_2,unit_3, unit_4, unit_5])); 
    
   cp_1100 = p1_1100/p_1100;   
    
p_1000 = mvncdf([-Inf, rf_cut_full(unit_2), rf_cut_full(unit_3),rf_cut_full(unit_4)], [rf_cut_full([unit_1]), Inf, Inf, Inf ],0,SIGMA([unit_1, unit_2, unit_3, unit_4],[unit_1,unit_2,unit_3, unit_4])); 
p1_1000 = mvncdf([-Inf,rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4),-Inf], [rf_cut_full([unit_1]), Inf, Inf, Inf, rf_cut_full(unit_5)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5],[unit_1,unit_2,unit_3, unit_4, unit_5])); 
    
   cp_1000 = p1_1000/p_1000;

p_0000 = mvncdf([rf_cut_full([unit_1]), rf_cut_full(unit_2), rf_cut_full(unit_3),rf_cut_full(unit_4)], [Inf, Inf, Inf, Inf ],0,SIGMA([unit_1, unit_2, unit_3, unit_4],[unit_1,unit_2,unit_3, unit_4])); 
p1_0000 = mvncdf([rf_cut_full([unit_1]),rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4),-Inf], [Inf, Inf, Inf, Inf, rf_cut_full(unit_5)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5],[unit_1,unit_2,unit_3, unit_4, unit_5])); 
    
   cp_0000 = p1_0000/p_0000;

cp_diff(1,4) = cp_1111-cp_0000;
cp_diff(2,4) = cp_1110-cp_0000;
cp_diff(3,4) = cp_1100-cp_0000;
cp_diff(4,4) = cp_1000-cp_0000;

   
%5th period 

p_11111 = p1_1111; 
p1_11111 = mvncdf([-Inf,-Inf,-Inf,-Inf,-Inf,-Inf], [rf_cut_full([unit_1]),rf_cut_full(unit_2), rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),rf_cut_full(unit_6)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 

    cp_11111 = p1_11111/p_11111;

p_11110 = mvncdf([-Inf,-Inf,-Inf,-Inf,rf_cut_full(unit_5)], [rf_cut_full([unit_1]), rf_cut_full(unit_2), rf_cut_full(unit_3),rf_cut_full(unit_4),Inf ],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5],[unit_1,unit_2,unit_3, unit_4, unit_5])); 
p1_11110 = mvncdf([-Inf,-Inf,-Inf,-Inf,rf_cut_full(unit_5),-Inf], [rf_cut_full([unit_1]), rf_cut_full(unit_2), rf_cut_full(unit_3),rf_cut_full(unit_4),Inf,rf_cut_full(unit_6)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 

    cp_11110 = p1_11110/p_11110;

p_11100 = mvncdf([-Inf,-Inf,-Inf,rf_cut_full(unit_4),rf_cut_full(unit_5)], [rf_cut_full([unit_1]), rf_cut_full(unit_2), rf_cut_full(unit_3),Inf,Inf ],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5],[unit_1,unit_2,unit_3, unit_4, unit_5])); 
p1_11100 = mvncdf([-Inf,-Inf,-Inf,rf_cut_full(unit_4),rf_cut_full(unit_5),-Inf], [rf_cut_full([unit_1]), rf_cut_full(unit_2), rf_cut_full(unit_3),Inf,Inf,rf_cut_full(unit_6)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 

    cp_11100 = p1_11100/p_11100;

p_11000 = mvncdf([-Inf,-Inf,rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5)], [rf_cut_full([unit_1]), rf_cut_full(unit_2),Inf,Inf,Inf ],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5],[unit_1,unit_2,unit_3, unit_4, unit_5])); 
p1_11000 = mvncdf([-Inf,-Inf,rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),-Inf], [rf_cut_full([unit_1]), rf_cut_full(unit_2),Inf,Inf,Inf,rf_cut_full(unit_6)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 

    cp_11000 = p1_11000/p_11000;
    
p_10000 = mvncdf([-Inf,rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5)], [rf_cut_full([unit_1]),Inf,Inf,Inf,Inf],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5],[unit_1,unit_2,unit_3, unit_4, unit_5])); 
p1_10000 = mvncdf([-Inf,rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),-Inf], [rf_cut_full([unit_1]),Inf,Inf,Inf,Inf,rf_cut_full(unit_6)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 

    cp_10000 = p1_10000/p_10000;
    
p_00000 = mvncdf([rf_cut_full([unit_1]),rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5)], [Inf,Inf,Inf,Inf,Inf ],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5],[unit_1,unit_2,unit_3, unit_4, unit_5])); 
p1_00000 = mvncdf([rf_cut_full(unit_1),rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),-Inf], [Inf,Inf,Inf,Inf,Inf,rf_cut_full(unit_6)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 

    cp_00000 = p1_00000/p_00000; 
    
cp_diff(1,5) = cp_11111-cp_00000;
cp_diff(2,5) = cp_11110-cp_00000;
cp_diff(3,5) = cp_11100-cp_00000;
cp_diff(4,5) = cp_11000-cp_00000;
cp_diff(5,5) = cp_10000-cp_00000;

%6th period 

p_111111 = p1_11111; 
p1_111111 = mvncdf([-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf], [rf_cut_full([unit_1]),rf_cut_full(unit_2), rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),rf_cut_full(unit_6),rf_cut_full(unit_7)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6, unit_7],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6, unit_7])); 

    cp_111111 = p1_111111/p_111111;
    
p_111110 = mvncdf([-Inf,-Inf,-Inf,-Inf,-Inf,rf_cut_full(unit_6)], [rf_cut_full([unit_1]),rf_cut_full(unit_2), rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),Inf],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 
p1_111110 = mvncdf([-Inf,-Inf,-Inf,-Inf,-Inf,rf_cut_full(unit_6),-Inf], [rf_cut_full([unit_1]),rf_cut_full(unit_2), rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),Inf,rf_cut_full(unit_7)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6, unit_7],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6, unit_7])); 

    cp_111110 = p1_111110/p_111110;
    
p_111100 = mvncdf([-Inf,-Inf,-Inf,-Inf,rf_cut_full(unit_5),rf_cut_full(unit_6)], [rf_cut_full([unit_1]),rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4),Inf,Inf],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 
p1_111100 = mvncdf([-Inf,-Inf,-Inf,-Inf,rf_cut_full(unit_5),rf_cut_full(unit_6),-Inf], [rf_cut_full([unit_1]),rf_cut_full(unit_2), rf_cut_full(unit_3),rf_cut_full(unit_4),Inf,Inf,rf_cut_full(unit_7)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6, unit_7],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6, unit_7])); 

    cp_111100 = p1_111100/p_111100;
    
p_111000 = mvncdf([-Inf,-Inf,-Inf,rf_cut_full(unit_4),rf_cut_full(unit_5),rf_cut_full(unit_6)], [rf_cut_full([unit_1]),rf_cut_full(unit_2),rf_cut_full(unit_3),Inf,Inf,Inf],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 
p1_111000 = mvncdf([-Inf,-Inf,-Inf,rf_cut_full(unit_4),rf_cut_full(unit_5),rf_cut_full(unit_6),-Inf], [rf_cut_full([unit_1]),rf_cut_full(unit_2), rf_cut_full(unit_3),Inf,Inf,Inf,rf_cut_full(unit_7)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6, unit_7],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6, unit_7])); 

    cp_111000 = p1_111000/p_111000;
        
p_110000 = mvncdf([-Inf,-Inf,rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),rf_cut_full(unit_6)], [rf_cut_full([unit_1]),rf_cut_full(unit_2),Inf,Inf,Inf,Inf],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 
p1_110000 = mvncdf([-Inf,-Inf,rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),rf_cut_full(unit_6),-Inf], [rf_cut_full([unit_1]),rf_cut_full(unit_2),Inf,Inf,Inf,Inf,rf_cut_full(unit_7)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6, unit_7],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6, unit_7])); 

    cp_110000 = p1_110000/p_110000;

p_100000 = mvncdf([-Inf,rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),rf_cut_full(unit_6)], [rf_cut_full([unit_1]),Inf,Inf,Inf,Inf,Inf],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 
p1_100000 = mvncdf([-Inf,rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),rf_cut_full(unit_6),-Inf], [rf_cut_full([unit_1]),Inf,Inf,Inf,Inf,Inf,rf_cut_full(unit_7)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6, unit_7],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6, unit_7])); 

    cp_100000 = p1_100000/p_100000;
    
p_000000 = mvncdf([rf_cut_full(unit_1),rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),rf_cut_full(unit_6)], [Inf,Inf,Inf,Inf,Inf,Inf],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6])); 
p1_000000 = mvncdf([rf_cut_full(unit_1),rf_cut_full(unit_2),rf_cut_full(unit_3),rf_cut_full(unit_4),rf_cut_full(unit_5),rf_cut_full(unit_6),-Inf], [Inf,Inf,Inf,Inf,Inf,Inf,rf_cut_full(unit_7)],0,SIGMA([unit_1, unit_2, unit_3, unit_4, unit_5, unit_6, unit_7],[unit_1,unit_2,unit_3, unit_4, unit_5, unit_6, unit_7])); 

    cp_000000 = p1_000000/p_000000;

cp_diff(1,6) = cp_111111-cp_000000;
cp_diff(2,6) = cp_111110-cp_000000;
cp_diff(3,6) = cp_111100-cp_000000;
cp_diff(4,6) = cp_111000-cp_000000;
cp_diff(5,6) = cp_110000-cp_000000;
cp_diff(6,6) = cp_100000-cp_000000;
